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ABSTRACT 

CN ; 

K*" , We formulate the collisionless Boltzmann equation (CBE) for dense star clus- 

00 . 

00 ■ ters that he within the radius of influence of a massive black hole in galactic 

■ nuclei. Our approach to these nearly Keplerian systems follows that of Sridhar 

■ and Touma (1999): Delaunay canonical variables are used to describe stellar 

orbits and we average over the fast Keplerian orbital phases. The stellar distri- 
bution function (DF) evolves on the longer time scale of precessional motions, 

. whose dynamics is governed by a Hamiltonian, given by the orbit-averaged 

^ ' self-gravitational potential of the cluster. We specialize to razor-thin, planar 
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discs and consider two counter-rotating ("±") populations of stars. To de- 
scribe discs of small eccentricities, we expand the ± Hamiltonian to fourth 
order in the eccentricities, with coefficients that depend self-consistently on 
the ± DFs. We construct approximate ± dynamical invariants and use Jeans' 
theorem to construct time-dependent ± DFs, which are completely described 
by their centroid coordinates and shape matrices. When the centroid eccen- 
tricities arc larger than the dispersion in eccentricities, the ± centroids obey a 
set of 4 autonomous equations ordinary differential equations. We show that 
these can be cast as a two-degree of freedom Hamiltonian system which is 
nonlinear, yet integrable. We study the linear instability of initially circular 
discs and derive a criterion for the counter-rotating instability. We then ex- 
plore the rich nonlinear dynamics of counter-rotating discs, with focus on the 
variety of steadily precessing eccentric configurations that arc allowed. The 
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stability and properties of these configurations are studied as functions of 
parameters such as the disc mass ratios and angular momentum. 

Key words: instabilities — stellar dynamics — celestial mechanics — galax- 
ies: nuclei 



1 INTRODUCTION 

Galactic nuclei have massive black holes and dense clusters of stars, whose structural and 



kinematic properties appear to be correlated wi ; 



1996 



Ferrarese fc Merritt 



2000 



Gebhardt et al. 



1 glob al galaxy properties (iGebhardt et al. 



20001) ■ These correlations are probably the 



relics of the formation and evolut ion of the galaxy and its central black hole ( iRichstone et al. 



1998 



Hopkins fc Quataertll201ll ). The dynamics of star clusters around massive black holes 



involves physical processes under ext reme conditions, because of th e high stellar densities, 



large velocities and short time scales (lAlexander 



2005 



Merritt 



20061 ). For a star cluster with 



1-dimensional velocity dispersion a, and black hole of mass M,, the radius of influence of the 
black hole is traditionally defined as rh = GM./o"^ . Within the radius of influence, i.e. for 
r < Th, the dynamics of stars is dominated by the gravitational attraction of the black hole. 
When general relativis tic effects are wea k enou gh, the orbital dynamics is a perturbation 



of the Kepler problem. 



Sridhar fc Toumal (119991 ) argued that the semi-major axis of stellar 



orbits would be a secularly conserved quantity, and that this greater integrability would 
facilitate the existence of asymmetric stellar distributions. The secular dynamical evolution 
of nearly Keplerian systems such as stellar clusters surr ounding black holes i n galactic nuclei, 



cometary clouds or planetesimal discs were studied by 



Touma et al. 



fl2009h . 



For most galaxies, it is difficult to observe details of the stellar distribution for r < 
Th ll Therefore, observations of the nuclear regions of our Galaxy and M31 assume special 
importance, because these are the nearest large, normal galaxies for which it is possible to 
get a great deal of photometric, kinematic and spectral information about the stars. There is 
evidence for a black hole of mass ~ 4 x 10® Mp^ at the Galactic center, with dense clusters of 



stars orbiting it ( 



Genzel et al. 



stars (jPaumard et al. 



20101 ). Among these, there is a population of about 200 young 
20061). about half which probably belong to a rotating disc which is 



highly warped (ILevin fc Beloborodov 



2003 



Lu et al. 



2006 



20091 ) . The remaining young stars 



appear to be members of a counter-rotating population which is thicker than and inclined 



1 rh ~ lOpc for a ~ 200 kms"! and M. ~ 10^ Mq. 
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to the warped disc (jGenzel et al. 



2003 



Paumard et al. 



2006 



Bartko et al. 



20091 . boioh . The 



nucleus of M31 has a lopsided double-peaked distribution of stars orbiting a black hole 
of mass ~ 10® Mfr^, with the brighter peak off-centered , and t he fainter one centered clos e 



to the black ho le (ILight et al. 



1974 : 



Lauer et al. 



199 



1998 



Kormendy &: Bender 



199% . 



Tremaind (Il995l ) proposed that the off-centered peak marks the location in a stellar disc 



corresponding to the aligned apoapsides of several eccentric stellar orbits. Following this 



proposal, more detailed stellar dynamica 


models of the eccentric disc have been pro 


posed 


(Bacon et al. 


2001 


Salow & Statler 


2001; 


Sambhus & Sridhar 


2002: 


Peiris & Tremaine 


2OO3I 



Bender et al 



20051 ). It is a very interesting fact that for both galaxies, these extraordinary 



stellar dynamical structures have r < r^. 

An alternative and useful way of seeing the significance of the radius of influence is 
this: a self-gravitating star cluster that lies within r^ has mass M < M,. The orbit of 
a star may be thought of as a slowly evolving Keplerian ellipse of fixed semi-major axis 
(with the central mass at one focus), whose dynamics is governed by a Hamiltonian which 



is th e orbit-averaged gravitational potential due to the star cluster (ISridhar &: Touma 



I999I ). This slow secular evolution time scale is larger than the typical Keplerian orbital 



time by the 



Sridhar et al. 



arge f a ctor (M,/M) . Sl ow modes of Keplerian discs were first explored by 
( 1999! ): iTremaind ( 1200 ll ). Orbital dynamics has b een analyzed and classified 



tials 


(Merritt & VaUuri 


1999; 


Sambhus & Sridhar 


2000 


Poon & Merritt 


2001 


Merritt & Vasiliev 


2011 


. One purpose of classifying stellar orbits is to be able to construct stellar distribution 



functions (DFs) that can reproduce the photometry and kinematics of stars around mas- 
sive black holes in galactic nuclei. When the time scales under consideration are smaller 
than the relaxation times, the DP obeys the coUisionless Boltzmann equation (CBE); see 



Binney fc Tremaind (l2008l ). In this paper we begin with a formulation of the CBE for the 



self-consistent, slow secular dynamics of star clusters within the radius of influence of mas- 
sive black holes. 

The warped di scs at the Galactic c e nter a re mutually counter-rotating. The stellar dy- 



namical model of 



Sambhus fc Sridharl (120021 ) for the lopsided stellar disc in the nucleus 



of M31 included a few percent of the stars on counter-rotating orbits. Here, it was pro- 
posed that the lopsidedness of the nuclear disc of M31 could have been excited by the 
counter-rotating instability, due to the accretion of a globular clust er that spiraled in due 
to dynamical friction. This proposal was motivated by the work of iToumal (120021 ). which 



Touma & Sridhar 



suggested that even a s mall fraction o f mas s in counter-rotating orbits could excite a lin- 



ear lopsided instability. 



Touma et al. 



( 120091 ) examine a secularly unstable system of coun- 



terrotating discs, and follow the unfolding and saturation of the instability into a global, 



uniformly precessing, lopsided (m = 1) mode. Counter-rotat ing streams o 



■ matter in a self- 



gravitating disc 


are known to be uns 


able to lopsided modes ( 


Zane & Hohl 


1978; 


Araki 


1987; 


Sawamura 


1988 




Merritt & Stiavelli 


199ol 


^almer & Papaloizou 


199ol: Isellwood & Merritt 


1994; 


Lovelace et al. 


1997; 


Sridhar & Saini 


2010 


). Mass accretion in galactic nuclei will be 



such that the sense of rotation of infalling material will be uncorrelated with the rotation 
of pre-existing material surrounding the central black hole. In other words, in the course of 
the evolution of a galaxy, having counter-rotating systems in its nucleus is probably generic. 

The main goal of this paper is to formulate and analyze the time-dependent dynamics 
of counter-rotating stellar discs, which lie within the radius of influence of the black hole. 
In § 2 we discuss the slow dynamics of nearly Keplerian star clusters by using the Delau- 
nay action-angle variables, and average over the fast orbital phase. We present the CBE 
governing the coUisionless evolution of the stellar DF. The self-consistent Hamiltonian is 
the orbit-averaged gravitational potential due to the star cluster, where softened gravity 
is used. The CBE is then formulated for razor-thin, planar discs, where the phase space 
is seen to be topologically equivalent to a 2-sphere. In § 3 we consider counter-rotating 
("±") discs with fixed semi-major axes. Here, it proves convenient to write separate CBEs 
for the ± DFs in the ± phase spaces. We consider discs of small eccentricities, and intro- 
duce "cartesian-type" canonical variables. The ring-ring interaction potential is expanded 



to 4t 
from 



1 order in the ecc e ntrici ties for discs with the same semi-major axes, using results 



Mroueh fc Toumal ( 1201 ll ) which are elaborated upon in Appendix A. Using this, the 



± Hamiltonians are expressed in terms of the ± DFs with coefficients that depend on both 
± DFs. In § 4 we construct time-dependent DFs for the ± discs. The method used is to 
first seek approximate dynamical invariants for the time-dependent dynamics of the gravi- 
tationally coupled ± discs, and then use Jeans' theorem to construct time-dependent DFs. 
The isocontours of the DFs are ellipses centered on moving origins in the ± phase spaces. 
The coordinates of the centers (referred to as "centroids") contain information about the 
mean eccentricities and periapse orientations of the ± DFs. Information about the disper- 
sions of eccentricities and periapse orientations is contained in the positive-definite 2x2 
"shape matrices" which describe the elliptical isocontours. When the centroid eccentricities 
are larger than the dispersion in eccentricities, the centroid dynamics is independent of the 
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shape dynamics; however, the shape dynamics is driven by the centroid dynamics. In § 5 we 
show that the coupled dynamics of the ± centroids is a two-degree of freedom Hamiltonian 
system which is nonhnear, yet integrable. We then study the hnear instabihty of initially 
circular discs and derive a criterion for the instability. The rest of the section is an explo- 
ration of the rich nonlinear dynamics of counter-rotating discs, with focus on the variety of 
steadily precessing eccentric configurations that are allowed. The stability and properties of 
these configurations are studied as functions of parameters such as the disc mass ratios and 
angular momentum. Summary and conclusions are offered in § 6. 



2 COLLISIONLESS EVOLUTION OF NEARLY KEPLERIAN STELLAR 
SYSTEMS 

We consider a stellar system around a central object of mass M, . Over times shorter than 
the relaxation times, the system is effectively collisionless. Then the stellar system may be 
thought of as composed of an infinite number of stars, each of infinitesimal mass, with total 
mass in stars equal to M. Stellar orbits are governed by the Newtonian gravity of the central 
mass, as well as the mean-field gravitational potential of all the stars. When (M/M,) ^ 1, 
it is useful to regard the dynamics as a perturbation of the Kepler problem. Thus the orbit 
of a star may be thought of as a slowly evolving Keplerian ellipse of fixed semi-major axis, 
with the central mass at one focus. This slow secular evolution time scale is larger than the 
typical Keplerian orbital time by the large factor (M,/M) . 

Each star is a moving point-like object in the 6-dimensional phase space, {r,v), where 
r is its position with respect to the central mass, and v is its velocity. Since the dynam- 
ics of a star is nearly Keplerian, it is useful to employ the Delaunay variables, which are 
action-angle variables for the Kepler problem. The Delaunay vari ables, |/, L,L^; w , q,h} 



are a set of action and angle variables for the Kepler problem ( ISridhar fc Touma 



1999; 



Binney fc Tremainell2008l ). The three actions are: / = \/GM,a where a is the semi-major 
axis; L, which is the magnitude of the orbital angular momentum; and Lz, which is the 
z-component of the orbital angular momentum. The angles conjugate to them are, respec- 
tively: w, the orbital phase; g, the angle to periapse from the ascending node; and h, the 
longitude of the ascending node. 

In the absence of self-gravity, the motion of the star is purely Keplerian: the orbital phase 
w advances steadily at a rate equal to the Keplerian orbital frequency, whereas the other five 
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Delaunay variables are constant in time. However, self-gravity contributes to a slow variation 
of the Delaunay variables. Let H be the total gravitational potential seen by a star, averaged 
over the Keplerian orbital phase of the concerned star. Then H = H{I, L, Lz, g, h, t), where 
we have allowed for a slow time dependence. Since H is — by definition — independent of 
w, the conjugate momentum, /, is a conserved quantity; the star's orbit can be imagined 
to be a slowly deforming "Gaussian ring" of fixed semi-major axis, with the central mass 
stationary at one focuso The slow secular evolution of the other Delaunay variables is given 
by 



dL 
d^ 



dH 

dg 



dg 
dt 



dH 
'dL 



dL, 



dH 
'dh 



dh 
dt 



dH 



[I) 



dt dh ' dt dLz 

Let the stellar system be described by a distribution function (DF), f{a,L,Lz,g,h,t), 
where / da dL dL^ dg dh is the mass in the element (da dL dL^ dg dh). The collisionless Boltz- 
mann equation (CBE) which describes the time evolution of the DF is 



where 



dt 



dl 
dt 



+ 



(2) 



df dH 
dg dL 



df dH df dH 
dL dg dh dLz 



df dH 



(3) 



dLz dh 

is the Poisson Bracket between / and defined in the (L, L^, (7, h) phase space. The Hamil 
tonian is 



H = <!>,,,ia,L,Lz,g,h,t) - G j da' dL' dL'^dg' dh' f{a',L',L'^,g',h',t) x 

X ^(a, L, Lz, g, h; a, L', L'^, g', h'), (4) 

where $ext is the gravitational potential due to an external source, averaged over the star's 
Keplerian orbital phase. The second term is the self-gravitational potential of the star clus- 
ter; the quantity. 



^ This is a restatement of the well-known result in planetary dynamics that the semi-major axis, a, is a secular invariant. 
Henceforth we use a instead of I. 
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^(a, L, L^, g, h; a', L', L'^,g', h') 



dw dw' 



1 



I 1/2 



(5) 



27r 27r j" |j. _ j./|2 _|_ ^2 j 
is proportional to the mutual gravitational potential energy between two particles, each of 
unit mass, averaged over their Keplerian orbital phases. Note that the gravitational interac- 
tion between the stars has a softening length b, whereas that between the central mass and 
a star i s thro ugh an unsoftened Keplerian potential. "Softened gravity" was introduced by 



Millerl ( 1l97ll ). and used as a surrogate for velocity dispersion in a stellar di sc; more recen t 



work using softened gravity in the context of nearly Keplerian systems are iToumal (|2002[ ): 



Touma et al. 



(120091 ). Since each star is usefully imagined to be a "Gaussian ring" in sec- 
ular dynamics, we refer to \E' as the ring-ring interaction function. At each value of the 
semi-major axis, a, equations ([2]) — (jl]) provide a self-consistent description of slow secular 
dynamics in the (L, Lz, g, h) phase space. 



2.1 The CBE for razor-thin discs 

When motion is restricted to a plane, the description of secular dynamics simplifies by 
reduction of the phase space by two dimensions. Let this plane be chosen perpendicular to 
the 2;-axis. The angles, g and h, no longer have clear independent meanings; rather it is 
the sum {g + h) that is well-defined, and we will henceforth refer to this as the angle g. 
Its conjugate variable is the scaled action variable, £ = L^j \/GM,a , which is equal to the 
2;-component of the angular momentum divided by the angular momentum of a circular 
orbit of radius a : this definition makes I a normalized quantity: — 1 ^ £ ^ 1 . Then 

dg_ _ dH_ cM _ _dH_ 

dt ~ ~dl' dt ~ ~~di' ^ ^ 

where the Hamiltonian H{a,i, g,t) = H/y/GM,a . The DF, f{a, i, g,t), satisfies the CBE: 



%.%.if,m.o, (7) 

where 

f OH f OH 
^^'^^ = 0^^ ~ Oi^ 
is the Poisson Bracket between / and H, defined in the {i, g) phase space. H is determined 

self-consistently through 
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H = $ext(a,£,^7,t) - ^J§^j da'dUg' f{a'A\g',t)-^{aA.9:a'A\9'). (9) 

Note that the kinetic energy and the Kepler part of the potential energy have been omitted 
because their sum which is equal to —GM,/2a is a constant. The {(!-■, g) phase space is, 
topologically speaking, a 2-dimensional sphere, with ^ equal to the cosine of the colatitude 
and g equal to the azimuthal angle: £ = ±1 are located at the north and south poles 
respectively, whereas £ = corresponds to the equator. For each value of the semi-major 
axis, a, equations ^ — ([9]) provide a self-consistent description of slow secular dynamics on 
this 2-sphere. 



3 COUNTER-ROTATING DISCS: FORMALISM FOR SMALL 
ECCENTRICITIES 

In our study of counter-rotating discs we restrict attention to isolated ($ext = 0) planar, 
razor-thin discs around a central mass. We also assume that all stars in a disc have the 
same semi-major axis: the + disc has stars with semi-major axis equal to a+ and — disc 
has stars with semi-major axis equal to a_. Thus we choose the DF to be of the form. 



f{a,e,g,t) = 6{a-a+)f+{e,g,t) + 5(a - a_) t), (10) 

where the (5-functions fix the semi-major axes of the two populations, which are now de- 
scribed by two different DFs, and /~. These DFs obey separate CBEs: 

^ + [r^H-] =0; ^ + [f-,H-] = 0, (11) 

where are the Hamiltonians acting on the ± populations, respectively. Each of and 
H depends on both /+ and /~, leading to coupled dynamics of the ± populations. In fact, 



H+{i,g,t) 



G 

^GM,a^ 



di'dg'f-^{i',g',t)<i/{a+J,g;a+,i',g' 



G 

^GM,a^ 



di'dg'f-{i\g',t)<il{a+J,g;a.,i',g' 



(12) 



and 
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H-{i,g,t) = - JL_ I dedg'f-{i',g',t)^ia^,i,g;a,J',g') 

G 



^^^rj^ j /+(^', 9\ t) vl>(a_, £, g- a+, f , g') . (13) 

We now specialize to counter-rotating discs of small eccentricities: let be the DF 
for a prograde population, which is concentrated around i = +1, and /~ be the DF for 
a retrograde DF which is concentrated around i = —1. We are interested in recovering a 
truncated model for the collective dynamics of these coupled populations. We recall that the 
{£, g) phase space is a 2-sphere, with I equal to the cosine of the colatitude and g equal to 
the azimuthal angle. The prograde and retrograde populations we consider are concentrated 
at the north and south poles, respectively. So it is convenient to use two different coordinate 
patches to describe the two populationsjj Thus we choose separate prograde and retrograde 



canonical variables 

/+ = e+ = -g- 

/_ = l + e_=g. (14) 

We will also find it convenient to use the "cartesian counterparts" of the (/, 9) variables. 
These are defined by 



x+ 



'2U sin^. 



-xnii 



sing . 



cosg; 



(15) 



X- = ^/2IZsm9- = ^2(1 + £) sing, 

y- = y2Jlcos^_ = ^/2{l + i) cosg. (16) 

Here x± are new coordinates, and y± are new momenta for the ± populations. The transfor- 
mations from old to new variable are of course canonical and can be simply recovered with 
the help of the generating function S{x±, 6±) = {x^/2) cot 9^-. 

^ Our analysis is limited to scenarios in which the two populations preserve their identity as prograde and retrograde stars. In 
other words, the sign of the orbital angular momentum of each star does not change. 
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Before we plunge into a series of approximations that will yield the reduced dynamics, 
we further simplify the model disc, by further specializing to ± populations with the same 
semi-major axes: 



a+ = a_ = ao . (17) 

The principal advantage of this restriction is that the ring-ring interaction function, can 
be described by fewer constants, allowing us to develop the theory with less clutter. However 
it may miss describing new phenomena when a+ 7^ a_. We reiterate: 

• We will study the planar secular dynamics of two counter-rotating stellar discs of small 
eccentricities, around a central massive object. 

• All the stars are assumed to have the same (conserved) semi-major axis. 

• Each star has a single degree of a freedom, namely the freedom to adjust its periapse ori- 
entation and conjugate eccentricity in response to collective (smooth) gravitational potential 
of all the other stars. 

• A star can, of course, change the sign of its orbital angular momentum and switch 
membership from the + population to the — population (or vice versa), while preserving its 
semi-major axis. However, this freedom is not allowed in what follows and shall be relaxed 
in future considerations of this problem. 



3.1 Expansion of the ring— ring interaction function 

To work out the Hamiltonians and we expand the ring-ring interaction func- 

tion, \1^, to 4*^^ order in the eccentricities of the ringcl. This expansion was developed by 



Mroueh &: Toumal (]201ll ). and is given in the Appendix A. Let us define e and e', the ec- 



centricity vectors characterizing two rings: 



e = (e cos (7 , e sin^f ) ; e' = (e' cos , e' sin ^f' ) , (18) 

with e = Vl — £2 and e' = \/ 1 — f I In the expansion of \I' we drop terms of the following 
type: (i) terms that are independent of e because these do not contribute to the dynamics 



^ Here a ring is thought of as a single Keplerian orbit which has been averaged over its fast orbital phase. In the context of 
stellar dynamics, we can also imagine a ring as a single Kepler orbit populated by many stars with the same Keplerian orbital 
elements, except for their orbital phases which are equally distributed over 2-k. 
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of the concerned ring; (ii) terms higher that 4*^ order in e, e', because this is the accuracy 
we aim for. Then, 



W— = ae2 + /3e-e' + 7e2e'2 + A(e-e')' + /s:(e-e')e'2 + xe' + «:e2(e-e') , (19) 
V M,ao 

where the coefficients (a, /3, 7, A, k, %) are functions of and b, and are given in terms 
of the softened Laplace coefficients (see the Appendix A). Note that 3'^'^ order terms are 
absent. Each of the e and e' can belong to either the + or — population, so there are 
four possibilities. We first express e in terms of {x±,y±) accurate to 4*^ order. For the + 
population equations (fT4l) and ( |T5l) give: 



= 1-^' = 1 - (1-/+)' = xl + yl- l{xl+yiy 



ecosg = W2/+ (^1 - COS0+ = y+y^l - i (a;^ + 



I4 



esmg = - ^/ 2/+ (^1 - sin = _ x+y^l - 1 (x^ ■ (20) 
Similarly, for the — population, equations (fT^ and (fT6|) give: 



ecos^ = J2I_ (^1 - cos9_ = y_^l - l(xl + y^J 



esing = ^/2/_ (^1 - y^ sin^^_ = x^^l - \{xl + yl) . (21) 

To express e' in terms of {x'^,y'^) we simply add the primes to the above expressions. Our 
next task is to obtain expressions for \l/(+,+'), \l/(+,— '), \l/(— ,+'), and \l/(— , — ') to the 
same orderjfl We begin with \l/(+,+'), by working out the terms involved in the (+,+') 
interactions. Using equation (pOj) . we have: 



^ We employ an obvious shorthand: for instance, '!'(+, — ') is the interaction function between two rings, one with parameters 
{x^,y-^.) and the other with parameters . 
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e = X 



+ y+ - -:{x+ + y+) 



2 , „,2\2 



e • e 



2 /2 

e e 



,/n2 



e ■ e 



e^e-e' 



+ y+y'+) \ (l 



x+x'^ + y+y'+) 



x\ + y\ 



x'l + y'l 



1 - 



xl + yl + x'l + y'l 



X 



+ yl) [x'l + y'l) + 



'^ + yW^T + ••• 



X 



-'^ + 1/'+) [x+x'^ + + 



X 



X 



2 + yl) {x+x'_^ + y+y'^) + 



(22) 



Then, substituting equations ( 12^ in (IT^ . the interaction function between two prograde 
rings reduces to: 



G 



+ 



+ 



+ 



+ 



+ 



a + 7 (^x'l + y'l) + \x' 
" + 7 [x'l + y'l) + Xy'l yl + [2As'^y '+] x+y^ 



i3x'+ + ( ~ ^) + + y 



/2 , /2 



x+ 



y+ 



K X _ 



/3\ , 
ti- -\y - 



x+ {xl + yl) 



y+{xl + yl) + (x-f) {xl + yl)' 



(23) 



where we have lumped all the dependences on primed quantities within the square brackets. 
Similarly, we compute quantities for the (+, — ') interactions using equations ( 12T|) . Then, 
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M.ao 



a + 7 ( x't + y''t ] + Ax'^ 



X, 



« + 7 [x'l + y'l) + \y't\ yl - [2\x'_y'_] x+y+ 



/3x'^ + (k — ^ j x'_ ^x'^ + y 



/2 , /2 



+ 



i3y'- + {^-^) y'- (^'- + y'- 



y+ 



/3 



K X _ 



a;+ + yl) 



+ 



3.2 Self-consistency: Hamiltonians in terms of the DFs 

We can now compute by using equations fl23|) and fl2^ in f|T2|) . Tlien, obtaining an 
expression for H is just a matter of switching signs: replace all the + variables by the 
— variables and vice versa. Putting together all these expansions, one finally recovers the 
Hamiltonians governing the prograde and retrograde populations: 



= \a^xI + 5±x±|/± + \c±yl + D±x± + E^y^ 



+ F^x^ {xl + yl) + [xl + yl)+K (x| + yl)\ (25) 



where the coefficients are determined self-consistently in terms of the DFs and / by, 
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A±it) = - 2 y" dx±dy± /±(x±, y±, [a + (7 + A)4 + 7y|] 

- 2 y dx^dy^ f^{x^, y^, t) [a + (7 + A)^ + 72/^ 

= -2Ay" da;±d2/±/^(a;±,2/±,t)a;±2/± + '^^ J dx^(^yTf^i^T^yT^'t)^TyT 

C±{t) = -2 J dx±dy± /^(x±, 2/±, t) [a + 7x| + (7 + \)yl] 

- 2 y dx^dy^ 2/^, [a + 7x| + (7 + A)2/|] 



-DiW = - / d2;±d2/±/=^(x±,2/±,^) 



/3a;± + (^'^ - + 2/|) 

+ J dx:^dy^f^{x^,y^,t) Px^ + (^k - X:^ {x^ + y^) 



E±{t) = - dx±dy±f^{x±,y±,t) 



-/ 



Py± + (^'^ ~ ^± + 



- K - 



where 



(x - f ) y dx+d2/+ r - {x-j)l dx_dy_ /- = - - ^) M , (26) 



± 



(27) 
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are the (constant) masses in ± populations, and M = M+ + M_ is the total mass in both 
discs. 

The coefficients [A^., B±, C±, D±, E^., F±, G±) are all, in general, functions of time, whereas 
i^' is a constant proportional to the total mass in both discs. Since we have defined separate 
canonical variables for the ± populations, it is necessary to take care to write the CBEs of 
equation (fTTI) as 

^ + [f^^HX = 0; ^ + [f-,H-]_ = 0, (28) 

where we have put ± subscripts on the Poisson Brackets to indicate that they are to be 
taken with respect to the appropriate set of canonical variables. Then the CBEs in equa- 
tion fl28|l together with the expressions for given in equations ( 125|) and fl26|) completely 
define the self-consistent evolution of the counter-rotating discs, accurate to 4*^ order in 
the eccentricities. 



4 COUNTER-ROTATING DISCS: TIME-DEPENDENT DFS 

In galactic dynamics, it is possible to construct many steady state solutions of the self- 
consistent CBE, whereas time-dependent behaviour is very difficult to understand even in 
the linearized limit. However, the present case of counter-rotating discs of small eccentric- 
ities around a central object turns out to be more tractable. The self-consistent dynamics 
described by equations (l25l) — (l28ll has implicit in it a certain approximate integrable dy- 
namics. This remarkable circumstance allows us to construct approximate time-dependent, 
self-consistent DFs, and describe the evolution of the counter-rotating instability largely 
analytically. 



4.1 An approximate dynamical invariant 

In this subsection we are interested in the construction of an approximate invariant for the 
dynamics on a two dimensional phase space, generated by a time-dependent Hamiltonian 
which is similar in form to those of equations ( l25l) . We drop the ± signifiers on all quantities 
in the interests of reducing clutter, but will restore them in the next subsection. Hence 
consider the Hamiltonian 
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H{x,y,t) = -A{t)x^ + B{t)xy + -C{t)y^ + D{t)x + E{t)y 

+ F{t)x {x" + y^) + G{t)y {x^ + y^) + K {x^ + y^Y . (29) 

This is a time-dependent Hamiltonian acting on the two-dimensional phase space {x,y). 
Wc now seek to eUminate the hnear terms in H by making a canonical transformation to 
new coordinate and momentum, {^1,^2), through a generating function 



Six,^2,t) = [x- + Y{t)] , (30) 

where X{t) and Y{t) are some time-dependent functions, which are to be determined. Since 
y — (dS/dx) and ^1 = {dS/d^2), the transformation 

X^i^ + X{t)- y = e2+n^), (31) 

amounts to a time-dependent shift of the origin of phase space. The new Hamiltonian is 
given by 

i/t(6,e2,i) = ^(^(6,e2),z/(a,e2),t) + ^ 

= i/un + + H['^ , (32) 

where i/iin, i/f^^ and H^^^ contain terms that are linear, quadratic, and cubic plus fourth 
order in (^1,^2), respectively. It is straightforward to work out that 



i/iin = AXi, + BYi^ + BXi2 + CYi2 + + E^^ 

+ Fei(x2 + r2) + 2FX(xei + r6) + G'6(^' + i^') + 2G'r(xei + f^s) 
+ [(r^ + rx2)6 + {X' + xr2)ei] + - ^6 ■ (33) 

We now require that i/iin vanishes. This happens when X{t) and Y{t) obey the following 
first order ordinary differential equations (ODEs): 
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dX 

= BX + CY + E + 2FXY + GiX"^ + SF^) + AKiY^ + YX^) 

at 

dY 

— = ~AX - BY - D - F(3X2 + f2) _ 2GXY - AK{X^ + XY^) . (34) 

Having eliminated i^un, the remaining terms in the new Hamiltonian are and h[^\ 
We require the former]^ 



^r(6,6,t) = + i?t(t)ei6 + \c,mh (35) 

where the new coefficients, 

— = - + 3FX + GY + 2ir(F2 + 3x2) 
Bt = B + 2FY + 2GX + 8KXY 

— = - + FX + 3GY + 2K{X^ + 3Y^) , (36) 
2 2 

are given in terms of the old coefficients and the centroid coordinates. 

The homogeneous, linear and time-dependent dynamics generated by preserves areas 
in (^i,.^2) space. Moreover, initial conditions given on any ellipse centered at = and 
^2 = will, at a later time, lie on some other centered ellipse of the same areajj Therefore, 
there must be a quadratic quantity that is preserved by the dynamics. Let us write this 
invariant as: 

■J = ifQit)^ = lQ^JitMJ, (37) 

where Q{t) is a time-dependent, positive definite, 2x2 matrix. Because phase areas are 
conserved, det(Q) is constant. The linear dynamics generated by H^^'' can be written in 
matrix form as: 



Tm; m = ^ (38) 

-A -B, 



dt 



^ The last bit of the new Hamiltonian, h[^^ = Ft£,i (Cf + €|) + Gt£,2 {£,1 +£.l)+K (^l + , contains cubic and fourth order 

terms in (^i , ^2 ) • The new coefficients, Ft = F + AKX and Gt = G + AKY , are given in terms of the old coefficients and the 

centroid coordinates. Henceforth we will not consider the modification of the dynamics due to these higher order terms. 
^ The shape — i.e. axis ratio — and orientation are determined by the matrix ODE Eg. 1501 
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where we have introduced T{t), which is a time-dependent, traceless, 2x2 matrix. If J, 
given in equation ( 1371) . is an invariant of the dynamics, we must have 



dJ ^ I t 
dt " 2^ 



T^Q + QT + 



dQ_ 



0, 



Therefore Q{t) must obey the matrix ODE: 

At 

It can be verified that the equation above preserves det(Q). 



T^Q - QT 



(39) 



(40) 



4.2 Distribution functions 

We are now ready to deal with the self-consistent dynamics of counter-rotating discs, de- 
scribed by equations (125!) — (!28|l . We now restore the ± signs that were dropped in the 
previous subsection. The first step is to pass from the Hamiltonians of equation f l25|) 
to new ± Hamiltonians h[^^^ which are of the form given by equation f l35p . We do not 
need to write down these new Hamiltonians; it suffices to note that they possess quadratic, 
time-dependent invariants of the form given in equations f l37|) : 



J±{x±,y±,t) = ]-QUt) [x± - X±{t)f + QUt) [x± - X±(t)] [y± - r±(t)] 



(41) 



The DFs f^{x±,y±,t) obey the CBEs of equations ( 1281) . By Jeans' theorem (IBinney fc Tremaine 



20081 ) . any function of the dynamical invariants is a solution of the CBE. Therefore we choose 



the ± DFs to be functions of the approximate invariants J± : 



f^{x±,y±,t) = F±(J±) . (42) 

Schematic representations of DFs in the {x±,y±) phase-spaces, and of centroid orbits in 
physical space, are shown in Fig.([T]). Here, we note the following important properties: 

(i) The DFs -F^(J±) are assumed to have compact support over the interval ^ J± ^ 

<-^itmax 1 • 

(ii) X±(t) and Y±{t) are the coordinates of the centroids of the DFs, {J±) , in the 
{x±,y±) phase spaces, respectively. 

(iii) At any instant of time, the isocontours of the DFs are ellipses in the {x±,y±) phase 
spaces that are centered on {X±{t),Y±{t)). 
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(a) y. 



(b) 





(c) 



(d) 




X_(() 



Figure 1. Distribution functions and centroid orbits: Panels "a" and "c" are schematic illustrations of the ± stellar distributions 
in the {x±,y±) phase spaces, respectively, at some given time. The elliptical patches indicate the regions populated by stars; 
each point in this phase space corresponds to a Keplerian orbit in physical space. The centroids of the ellipses, (X±{t), Y±(t)), 
are marked by black dots, and the principal axes of the ellipses are given by the eigenvectors of the symmetric matrices, (t) . 
Panels "b" and "d" are representations, in physical space, of the Keplerian orbits corresponding to the ± centroids. Note that 
the angles to the pericentre, = =p arctan(X±/Y±) ; hence, for the distributions shown in panels "a" and "c", we have < 
and gi > . 



(iv) The shapes and orientations of the elhpses are described by the time-dependent, 
symmetric, positive definite matrices Q^{t), which we will refer to as the shape matrices. 
Time evolution preserves det (Q±); we can choose det (Q^) = 1 ■ 

(v) The zeroth and first moments of the DFs, 



M4 



27idJ±F^{J±) ; 



a. 



± 



J±n 



(43) 



are the disc masses and the squared dispersions (of eccentricities), respectively. Once the 
DFs have been specified, M± and cr^ can be treated as constants. 

(vi) We can now state precisely the conditions under which -F^(J±) are approximate 
solutions. These DFs are good solutions when the dispersions in the eccentricities are much 
smaller than their centroid values; i.e. when 



^ o-± < e± < 1 



(44) 



where e± is a typical value of a/ X'^ + . 

(vii) The total angular momentum in the two discs is 

Aot = VGM.ao j didgif-^ + ^GM.ao j d£dg i T , 

where £ is positive for /+ with didg = dx+d?/+, and i is negative for /~ with d£dg = dx^dy^. 
Using equations ffT4l) — (fT6|) and (I3T]) . we write 



(45) 



± 



1 - ^(4 + y|) 



± 



{(6± + X^f + iU + Y^f} 



(46) 
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For the DFs given in equation P2|) . /^(x±, y±,t) = {J±)- When equation f l46|) is used, 
and ^2± contribute to the integral of equation (145|) only at second order. These contributions 
are small, since the dispersion in eccentricities is much smaller than centroid values for these 
DFs. So we drop the dependences on ^i-t and ^2± on the right side of equation ( l46ll . and 
set x± = X± and y± = Y± . Using the definitions of M± given in the first of equations (l43ll 
above, we obtain 



- M_] - ^[Xl + Yl] + ^ [X^ + Y'] + O {al/el) . (47) 



VGM,ao 



The first term on the right side is the contribution from the ± discs if they were circular; 
the second term is the decrement due to the centroid eccentricity of the + disc; the third 
term is a similar and oppositely signed contribution from the — disc. 



We need to compute the coefficients {A±, B±, C±, D±, E±, F±, G±, K± ), by substituting 
equation ( H2l) for the ± DFs in equations ( l26l) . Similar to the treatment of the angular mo- 
mentum of the discs given above, we set x± = X± and y± = Y^. . Then, it is straightforward 
to express the coefficients as functions of the centroid coordinates, [X±(t) , Y±(t)]: 
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A± ^ - 2Ma - 2M± [(7 + A)X| + ^Y^] - 2M^ [(7 + \)X^ + jY^] 



B± = -2\M±X±Y± + 2\M^X^Y^ 



C± = - 2Ma - 2M± [^Xl + (7 + A)r|] - 2M^ [^Xl + (7 + A)r|] 



- - M± + x±y|] + ^M^x^ + ("^ " + 



± 



- I « - ^ ) [M±X± - il^X^] 



[M±y± + M^Y^] 



The centroid coordinates obey: 



(48) 



= s±x± + c±y± + E± + 2F±x±y± + G± {xl + syj) + ak± (yj + y±x|) 

dt 



-A±X± - B±Y± -D±-F± {3X1 + y|) - 2G'±X±y± - 4K± (X| + X±yJ) 



(49) 

These are a set of 4 autonomous first order ODEs with cubic nonhnearity. The shape matrices 
obey the following first order matrix ODEs: 



dQ 

~dt 



± 



\-At± - Bt± ) 



(50) 



where 
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At± = A± + 6F±X± + 2G±F± + AK± (fJ + 3X|) 
Bt± = S± + 2F±y± + 2G±X± + 8/r±X±F± 

Ct± = C± + 2F±X± + 6G±y± + 4i^± (X^ + 3Y^) . (51) 

Since T'^{t) are traceless matrices, det [Q''^(t)] is a conserved quantity; without loss of 
generality, we can choose det [Q^(0)] = 1. The matrix T depends on the centroid co- 
ordinates, so the matrix equations for Q, while linear, are driven by centroid evolution. 
Equations (HHl) — fl5T]) determine the self-consistent centroid and shape dynamics of DFs 
describing the counter-rotating discs. 



5 COUNTER-ROTATING DISCS: CENTROID DYNAMICS 

As shown above, in the limit ^ cr± ^ e± ^ 1 , shape dynamics is driven by centroid 
dynamics, and consists of area preserving evolution of the shape and orientation of the 
elliptical isocontours of the DFs, with no feedback on centroidsjj In what follows, and with 
the understanding that shape can be recovered easily from centroids as and when required 
by a given application, we drop any further reference to shape dynamics, and focus our 
attention on the nonlinear evolution of the centroids. 



5.1 Integrability 

The centroid equations ( l49l) are a set of 4 autonomous first order ODEs with cubic nonlin- 
earity. Quite remarkably, it turns out that they describe a non linear, yet integrable, system. 
This happens because of the underlying Hamiltonian structure and the presence of a second 
conserved quantity. Let us define new rescaled variables (M±,f±): 



u± 



so + /i„) = 1 . 



(52) 



Such a feedback requires a higher order theory, and may ver y well account for instability saturation in a planar analog of the 



three dimensional saturation described in 



Touma et al. 



1 2009) 
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A lengthy yet straightforward calculation shows that the centroid equations fj49l) are equiv- 
alent to the following two degree-of-freedom system Hamiltonian system, with coordinates 
u± and momenta v± : 



dn 

dv^ ' 



di;+ 
'dt 



m 



dii- 
~dt 



m 

dv^ 



dv- 



m 



n 



+ Ma) [ul + vl] - + Ma) [ul + vl] 



7 + A 4x - a 
2 4/i+ 



2 4/i_ 



AM [u+u- - v+v-f - 7M [[u\ + (m?. + t;!)] 



+ ( K — — ) M [m+M- — w+f _] 



/ / 2 I 2 \ I / A'"- / 2 I 2 \ 



(53) 



where m_, f+, f+) is the Hamiltonian for the two degree-of-freedom system with 

coordinates u± and momenta v± . Since "H is independent of time, it is conserved. It is 
straightforward to verify that equations f l53|) also conserve the total angular momentum 
defined in equation ( |47|) . Dropping the first term, we write the second conserved quantity as 



C 



(54) 



4 4 

This quantity is a measure of the amount by which the angular momentum is lower than 
the maximum value it can attain when the centroid eccentricities are zero; for brevity we 
shall henceforth refer to C as the angular momentum. 

Since this two degree-of-freedom system has two independent conserved quantities, l-L 
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and £, it is integrable. We will explore the non linear dynamics of this system, after exam- 
ining the linear instability of zero eccentricity discs. 



5.2 Linear instability of zero eccentricity discs 

The zero eccentricity state, u+ = v± = Q which has £ = , is an equilibrium state of 
the dynamics governed by Ho Here we determine the conditions under which this equilib- 
rium is unstable to small perturbations. The linearized equations obeyed by infinitesimal 
perturbations are: 

dii+ dv^ 

— W+f + — WcV- , — - — = — WcU- ; 



dt ^ ^ ' dt 

du_ dv_ 

— — = — — W/U+ , — — = w^u_ — WcU+ , (55) 

dt ^ dt ^ ^ ^ 

where w± = {2Ma + (3M±) , and Wc = fiyJMj^M^ are constants. It is readily verified 
that this linearized system conserves C To solve these equations let us define the complex 
variables: 



2+ = + i f + , 2;_ = f _ + i M_ , (56) 

in terms of which equations f l55|) reduce to: 

dz+ . dz- . 

= 1W+Z+ - WcZ-, = -1W-Z_ - WcZ+. (57) 

Note that the asymmetry in Eqs. (1571) is inherited from the asymmetry in the definition of 
z± in Eqs. ( 156|) . Looking for normal modes, z± = Z±e^^ , we obtain and solve a quadratic 
characteristic equation for the two eigenvalues: 



s = 2^^^+"^-) ^ 2V^^^ ~ (w+ + w_) . (58) 
There is a growing solution when Aw"^ > (w_|_ + w_)^. Thus the zero eccentricity equilibrium, 
u± = v± = , is unstable (overstable) when 

^ Implicit to this zero eccentricity equilibrium are DFs of the form F^(J±) = ~^~^{J±): i-G, DFs with a"^ = and no shape 
(Q^) dynamics to speak of. 
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0.5 1 1.5 2 0.2 0.4 0.6 0.8 1 

r H 

Figure 2. Linear Stability of the zero angular momentum equilibrium: Left: Curve of critical values in the plane of mass ratio 
fi and dimensionless softening r = b/ao above which the centroids of the prograde and retrograde population are unstable to 
eccentricity growth. Right: The real parts of the eigenvalues given in equations I I58II . as a function of fi for r = 0.1 . Instability 
sets in for fi ^ ^tcrit — 0.4436 . 



Using the definitions of a and /3 given in tlie Appendix A, we display tlie stability criterion of 
equation (159|1 in Fig.® by plotting the mass ratio, /i = M_/M+ , versus the dimensionless 
softening, r = 6/aoo For a given softening, there is a critical value of the mass ratio 
(which decreases with increasing softening) above which the zero-eccentricity equilibrium 
is unstable. Conversely, for a given mass ratio, there is a critical value of the softening 
(which increases with decreasing mass ratio) above which the zero-eccentricity equilibrium 
is unstable. In the right panel of Fig.(|2]) we plot the real parts of the eigenvalues given in 
equations fl58|) . as a function of /i for r = 0.1 . When /i ^ yUcrit — 0.4436 , the eigenvalues are 
both imaginary, corresponding to normal modes describing steadily precessing discs of fixed 
centroid eccentricities. Both growing and damped solutions are allowed when fi > yUcrit • The 
growing (damped) solution describes discs whose eccentricities grow (damp) as they precess 
steadily. 

Since C is conserved, the instability operates through exchange of angular momentum 
between the prograde and retrograde discs. When the prograde disc gives some angular 
momeiitum to the retrograde disc, both discs increase their eccentricities. As is well-known 



flMiller 



1971 



Binney fc Tremaine 



20081 ). the softening length mimics the epicyclic radius of 



stars on nearly circular orbitso For this exchange to be self-reinforcing, for the instability 
to kick in, the mass ratio has to be large enough for ± disk self-gravity to overcome the 



Without loss of generality, we assume that ^ fi 1 . 

Our general formalism does not of course need softening. But, the DFs considered in Sec l4.2l and after are cold in the 
sense that the dispersions in the eccentricities are much smaller than the centroid eccentricities. So, softening serves to mimic 
eccentricity dispersion in our model (just like it mimics velocity dispersion in disc dynamics). 
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effective "heat" due to softening; ( a process which is similar to the one driving the radial 
orbit instability ( iLynden-Belll Il979h ) . 



5.3 Nonlinear dynamics 

As discussed earher centroid dynamics is integrable, because it is a two degree-of-freedom 
system (4-dimensional phase space) with two independent conserved quantities. We now 
use the conservation of the angular momentum of equation fl5^ to convert the problem 
into a Hamiltonian system of one degree-of-freedom system. We achieve this through two 
canonical transformations to new variables. First, we pass from {u±^ v±) to new action-angle 
variables, {L±^il)±): 



a/2L+ sin ijj^ 



\/2L+ cos ip^ 



M_ = a/2L_ sinip- , V- = ■\/2L^ costjj- ■ (60) 

Written in terms of the new variables, the Hamiltonian of equations f l53p becomes: 

H = — W-L^ — 2wc^/L+lI cos + ip^) 



■ 4:r]+Ll - 4t]_L^_ - 4 ^/t - M^L+L_ cos (V^+ + 





-4AML+L_cos2 (?/>+ + ?/'_) - 47ML+L_ 
where the new constants, 7]^, are defined by: 



(61) 



The angles ■?/'+ and ip- appear only in the combination {ipj^ + ipS). So we transform 
to new action-angle variables, (S, G) and {C,d), defined through the generating function, 
S (£, S, 7/^_) = (V^+ + 7/^„) S + {ip+ - ^_) £ . Then, 



•d = — ip- ^ Q = ip+ + ip- . 

When expressed in terms of the new variables, the Hamiltonian becomes: 



(63) 
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n 



(w+ + w^) Tj — {w^ — wJ) C — IWc^Y? — C? cos 9 



/3 



- 4 U - - MVS2 - £2 



cos B 



- 47M (E^ - _ 4^xM (S^ - £2) cos20 _ (64) 

Since H is independent of i), its conjugate momentum C is conservedo Hence C may be 
treated as a constant parameter occurring in H, which can be thought of as a Hamiltonian 
describing the one degree-of-freedom Hamiltonian with coordinate G and momentum E. 

The global structure of dynamics in the (O, S) phase space is most easily visualized by 
plotting the level curves of "H for different values of the constant parameter CilB It seems 
simplest to label the axes of the figures, using the "cartesian-type" canonical variables: 

U = V2(S - £) sine, V = ^2 (S - £) cosG, (65) 

where U and V are new coordinates and momenta, respectively. Before we discuss the phase 
space structure, it is useful to interpret the canonical variables in terms of physical quantities 
related to the centroids of the prograde and retrograde populations: 



• \/t/2 + V^2 _ 1^2(2 — C) = y/2LZ = + ^2 is proportional to the centroid eccen- 
tricity of the retrograde ring. 

• U = V = represents the zero eccentricity state; it is an equilibrium for £ = 0, but 
not otherwise. 

• O = {^|J-^- + ip^) = (^g^^^ — (7^''°*) is the difference between the periapse angles of the 
centroids of the retrograde and prograde populations. Since we also have 6 = arctan {U /V), 
we note that: 

(i) U = and V < implies that ((^f^"* — 51'^®°*) = tt, corresponding to eccentric ± discs 
with anti-aligned periapses. 



From equations II63I I and II6OI I. we have C = (L+ — L—) /2 = (n^ + ti^ — u'^ — v^) /4 equal to the conserved quantity 

defined earher in equation II 541 1. 

We have also confirmed that the same results are obtained through direct integration of the equations of motion. 
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(ii) f/ = and V > Q implies that (^g"^^^^ — (7^°"*) = 0, corresponding to eccentric ± discs 
with aligned periapses. 

Before we begin, we need to determine the parameters — other than C - that control the 
dynamics. We first note that the constants w±, Wc and ri± are all proportional to the total 
mass, M, in the ± discs. Each of the parameters (a, /3, 7, A, /t, x) is proportional to Oq^^^, 
where oq is the semi-major axis; these parameters are also functions of the dimensionless 
softening parameter r = b/aQ, but the dependences are not so simple. Therefore, every 
term on the right side of equation is proportional to Mqq ^^"^ , so this combination 
of constants determines only the rate at which a phase trajectory is traversed. For the 
purposes of investigating the geometry of phase space, we may set MaQ^^"^ equal to unity. 
The dimensionless masses, /i+ = {Mj^/M) and /i_ = (M^/M), can both be expressed in 
terms of the dimensionless mass ratio, fi = (M_/M+) . Therefore, we are left with the three 
dimensionless parameters {C,r,fi), which control the dynamics generated by "H . 

5.3.1 Dynamics when C = 

When £ = 0, the prograde and retrograde discs have equal amounts of mass-weighted 
centroid eccentricities. We have already studied the linear instability of the zero eccentricity 
equilibrium in § 15.21 At a fixed value of r, there is a critical value of the mass ratio fi above 
which the instability sets it. When r = 0.1, this critical mass ratio is ficnt — 0.4436, which 
corresponds to about 70% of the disc mass in the prograde component and the remaining 
mass in the retrograde component. Here we explore the structure of nonlinear centroid 
dynamics as fi is varied, with both C and r held fixed. 

Fig. ([3]) displays the level curves of "H in the {U, V) phase space for four different values of 
/i, when C = and r = 0.1 . Phase flows occur along these level curves. Some noteworthy 
features are: 

• For fi ^ /icrit, the zero-eccentricity equilibrium P2 = (0, 0) is stable; but there are also 
two additional equilibria. Pi = (0, Vi) and P3 = (0, V3). 

• Pi is stable: it has U = and Vi < 0, which corresponds to steadily precessing eccentric 
± discs with anti-aligned periapses. 

• P3 is unstable: it has U = and V3 > 0, which corresponds to steadily precessing 
eccentric ± discs with aligned periapses. 

• When /i exceeds yUcnt, the zero-eccentricity equilibrium P2 goes unstable by merging 
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Figure 4. Snapshot of the it centroid orbits in physical space corresponding to the stable, precessing equilibrium PI {U = 
0, V ~ —0.04) in the bottom-right panel of Fig. ((Sjl- The cross marks the location of the central mass, which is at the common 
focus of both the prograde (solid) and the retrograde (dashed) centroid orbits. The eccentricities of the ± orbits are 0.05 and 
0.06, respectively, and have been greatly exaggerated in the figures. Note that the ± periapses are anti— aligned. 

with the unstable equihbrium P3. Small perturbations about P2 now exhibit large variations 
in eccentricity, with phase space trajectories that take them on an excursion around the 
stable equilibrium Pi . A schematic view in physical space of anti-aligned centroid orbits at 
Pi is shown in Fig.(j4l). 

• Pi remains stable for all values of fi . 

• The bifurcation of equilibria as a function of /i, at fixed r, is shown in Fig. ([5]). The 
plot reflects transition across /icrit? with the two equilibria P2 and P3 merging to form an 
unstable equilibrium, along with the continuing sequence of the stable equilibrium Pi . 
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Figure 5. Bifurcation of equilibria: Plot of the values of the V-coordinates of the stable and unstable equilibria, as functions 
of the mass ratio fi. Stable/unstable equilibria are shown in solid/dashed lines. 

• Further increase in fi does not alter the equihbrium structure. However, there are larger 
excursions in eccentricity, so much so, that the conditions under which the model works can 
be in question. 

5.3.2 Dynamics when C ^ 

When £ 7^ 0, the mass-weighted centroid eccentricities of the prograde and retrograde discs 
are unequal. Without loss of generality, we assume that £ ^ 0; in other words, we assume 
that the mass-weighted centroid eccentricity of the prograde disc is greater than that of the 
retrograde disc. We follow the equilibria and their bifurcations as £ is increased, at fixed fj, 
and r. Some noteworthy features of the phase-space evolution, shown in Fig. ([6]), are: 

• As we have seen earlier, when C = , there is a qualitative change in the phase portrait 
when fi is smaller or larger than /icrit , for some chosen value of r. For r = 0.1,/i = 0.1< yUcrit 
so the zero-eccentricity equilibrium: we have three equilibria; two stable (Pi, P2) and one 
unstable (-P3), as discussed above. 

• With increasing C, P2 (which was initially at the origin) shifts continuously to higher 
eccentricity with U2 = and V2 > , corresponding to steadily precessing eccentric ± discs 
with aligned periapses. At the critical value of £ ~ 0.0271 , P2 and P3 collide. 

• All through. Pi remains stable with Ui = . However, as £ increases, Vi also increases 
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Figure 6. Phase-space evolution with increasing value of C. 



(a) (b) 




Figure 7. Snapshots of the ± centroid orbits in physical space, corresponding to the stable (precessing) equilibrium PI for 
two different cases. The cross marks the location of the central mass, which is at the common focus of both the prograde 
(solid) and the retrograde (dashed) centroid orbits. Panel "a" corresponds to top-right panel of Fig. (|6]l {C= 0.0272), when PI 
{U = 0,V ~ 0.07) is still a stable equilibrium before bifurcation. The eccentricities of the ± orbits are 0.35 and 0.2, respectively. 
Note that the it periapses are aligned. Panel "b" corresponds to the bottom-right panel of Figure 4 (£= 0.0272), when PI 
has bifurcated into two stable equilibria, Pla and Plb, which are at (7 ~ ±0.04 and V ~ 0.18; we show orbits for Pla . The 
eccentricities of the ± orbits are 0.65 and 0.55, respectively. Note that the it periapses are mis-aligned, with the pericentre of 
the "-" orbit leading the pericentre of the "+" orbit by about 12.5 deg (for Plb, the pericentre of the "-" orbit would lag the 
pericentre of the "+" orbit by the same amount. 

and, near C ^ 0.0031 , Vi becomes positive from its initially negative value; thus the cor- 
responding stable, steadily precessing eccentric ± discs switch from anti-aligned to aligned 
periapses. 

• With further increase in C, Pi remains stable with Ui = and increasing Vi until we 
hit another critical value, C ^ 0.0433 . At this value of C, Pi becomes unstable and, in a 
pitchfork-like bifurcation, there emerge two stable and non-aligned equilibria. The remark- 
able feature of these new equilibria is that they are neither aligned nor anti-aligned. This 
suggests that, for large enough values of £, the stable, uniformly precessing counter-rotating 
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Figure 8. Uniformly preccssing configurations for C . Left: Plot of the y-coordinatcs of the equilibria as functions of C, for 
the on-axis {U = 0) equilibria which correspond to aligned and anti-aligned periapses. Stable/unstable equilibria are shown in 
solid/dashed lines, unstable. Two critical values of C are apparent: C ~ 0.0271, at which stable and unstable equilibria merge; 
C ~ 0.0433, at which Pi becomes unstable, giving birth to non-aligned equilibria (with U ^ 0). Also apparent is the transition, 
at £ ~ 0.0031, of the Pi, from being an equilibrium with stable anti-aligned periapses (Vi < ) to stable aligned periapses 
(Vi > 0). Right: The bifurcation of the Pi, from being an aligned on-axis equilibrium to non-aligned equilibria. 

discs have periapses that are neither ahgned nor anti-ahgned. The transition from ahgned to 
non-ahgned stable equihbria at the bifurcation is illustrated in Fig.(I7j) with centroid orbits 
in physical space, at Pi before the bifurcation, and at one of its two stable offsprings after. 

• Increasing C still further maintains the equilibrium structure with the stable equilibria 
increasing in eccentricity and misalignment, in a fashion reminiscent of the displacement of 
an unstable bead on a rapidly spinning hoop. 

• In the left panel of Fig. ([8]), we show the on-axis equilibria, the eventual merging of P2 
and P3, followed by the loss of stability of Pi. In the right panel, we follow the bifurcation 
from Pi (now unstable) of the two stable, non-aligned equilibria. 

6 SUMMARY AND CONCLUSIONS 

We have formulated the problem of the coUisionless, self-consistent dynamics of nearly Kep- 
lerian star clusters around a massive black hole. Averaging over the fast Keplerian phase, we 
used the result that the semi-major axes of the stellar orbits are nearly conserved quantities. 
Hence each stellar orbit may be imagined to be a Gaussian ring, of fixed semi-major axis, 
whose shape and orientation changes over time scales that are longer than the Kepler or- 
bital times. Since the semi-major axis of each stellar orbit is constant, and the orbital phase 
unimportant, the phase space in which the slow dynamics of an individual orbit occurs is 
four-dimensional. In terms of the Delaunay actions-angle variables, the two actions are the 
magnitude of the angular momentum and the ^-component of the angular momentum, with 
corresponding angles being the argument of the periapse and the longitude of the ascending 
node. The star cluster was described by a distribution function (DP), which is a function of 
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five variables: the four action-angle variables described above, and the constant semi-major 
axis which can differ from star to star. We presented the collisionless Boltzmann equation 
(CBE), which governs the self-consistent slow dynamics of the star cluster. The Hamiltonian 
determining this dynamics is the orbit-averaged gravitational potential energy between the 
stars; it is determined by integrating a softened Keplerian potential over the orbital phases 
of both rings (thereby forming a "ring-ring interaction function"), and then integrating this 
over all of phase space, suitably weighted by the DF. 

The goal of this paper is to explore the counter-rotating instability. To this end, we 
considered the CBE for razor-thin, planar discs; in this case, at each value of the semi-major 
axis, the phase space is two dimensional, topologically equivalent to the 2-sphere. Then we 
considered the counter-rotating discs as two separate collisionless populations of stars, the 
prograde (or "+") disc and the retrograde (or "— " disc). For convenience, we assumed that 
all the prograde stars have the same semi-major axes and all the retrograde stars have the 
same semi-major axes. Thus the total phase space is the direct sum of two 2-spheres; one 
for the prograde disc, the other for the retrograde disc. We then wrote down two separate 
± DFs which obeyed two separate ± CBEs, governed by two separate ± Hamiltonians. 
Since the ± populations are gravitationally coupled, each of the ± Hamiltonians depends on 
both the ± DFs. When the discs are composed of stellar orbits of small eccentricities, the 
prograde population is clustered around the north pole of its 2-spherical phase space, and 
the retrograde population is clustered around the south pole of its 2-spherical phase space. 
We transformed to new "cartesian-type" canonical variables, which are more convenient to 
use in the limit of small eccentricities. The ring-ring inte raction function was then expanded 
to fourth order in the eccentricities, using results from iMroueh fc Toumal (120111 ). So, each 
of the ± Hamiltonians was obtained as a quartic polynomial of the ± canonical coordinates 
and momenta, with coefficients that depend on both the ± DFs. 

Time-dependent DFs were constructed using Jeans' theorem. The first step was to iden- 
tify appropriate phase space functions which are approximately conserved by the time- 
dependent ± Hamiltonians. Canonical transformations to moving origins in the ± phase 
spaces were used to construct these approximate dynamical invariants; to lowest order in 
the eccentricities, these are quadratic functions of the phase space coordinates and mo- 
menta with time dependent coefficients. The next step was to choose the ± DFs to be 
some physically allowed function (positive, finite mass, low eccentricity stars etc.) of the 
± invariants. The DFs are such that their isocontours in the ± phase spaces are ellipses. 
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centered on moving origins. The time-dependent coordinates of the origins were referred 
to as the centroids. The evolving shapes and orientations of the elhpses were described by 
± shape matrices, which are time-dependent and positive-definite; these matrices describe 
the anisotropic nature of the dispersion in the eccentricities. It is important to note that 
the DFs so constructed are approximate time-dependent solutions, when the dispersions of 
eccentricities described by the shape matrices are less than the centroid eccentricities. Thus 
the coupled, time-dependent dynamics of the the counter-rotating discs is described by a 
set of ordinary differential equations (ODEs) describing the evolution of the centroids and 
shape matrices. Some general properties are: 

(i) For each of the ± populations, there are two centroid coordinates and one positive- 
definite 2x2 shape matrix. Thus we have a system of 12 first order ODEs describing 
counter-rotating disc dynamics. 

(ii) The centroid equations are a set of 4 autonomous (i.e. time-independent) first order 
ODEs. These are independent of the shape matrices, because we are working in the limit 
where the centroid eccentricities are greater than the dispersion of eccentricities. 

(iii) The centroid equations conserve the total angular momentum corresponding to the 
centroid eccentricities. 

(iv) The equations for the shape matrices are linear ODEs which depend on the centroid 
coordinates. Thus shape dynamics is driven by centroid dynamics. 

(v) The shape equations conserve the determinants of the shape matrices, so we may 
choose the determinants to be equal to unity. 

(vi) The 4 autonomous first order ODEs describing centroid dynamics have cubic non- 
linearity. However, quite remarkably, they constitute an integrable system. This happens 
because the 4-dimensional system corresponds to a 2-degree of freedom Hamiltonian system 
which admits two conserved quantities; the centroid angular momentum and the Hamilto- 
nian itself. 

We then studied the linear stability of initially zero-eccentricity discs, and derived the 
conditions under which the configuration is unstable. Some notable properties are: 

(i) For a given softening, there is a critical value of the mass ratio (which decreases with 
increasing softening) above which the zero-eccentricity equilibrium is unstable to the growth 
of eccentricities in both ± discs. 

(ii) Conversely, for a given mass ratio, there is a critical value of the softening (which 
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increases with decreasing mass ratio) above which the zero-eccentricity equihbrium is un- 
stable. 

(iii) The stable solutions correspond to normal modes describing steadily precessing discs 
of fixed centroid eccentricities. 

(iv) When the parameters are in the unstable regime, both growing and damped solutions 
are allowed. The growing (damped) solution describes discs whose eccentricities grow (damp) 
as they process steadily. 

(v) The physical basis of the instability is through exchange of angular momentum be- 
tween the ± discs. When the prograde disc gives some angular momentum to the retrograde 
disc, both discs increase their eccentricities. For the instabihty to operate, the mass ratio 
has to be large enough to be able to overcome the effective "heat" in the stellar distribution 
which, in our cold DFs, is mimicked by the softening. 

The nonlinear dynamics is, of course, much richer. Wc demonstrated that the ODEs 
of centroid dynamics could be cast into Hamiltonian form, with a 2-degree of freedom, 
time-independent Hamiltonian which is quartic in the canonical variables. This Hamilto- 
nian conserves two independent phase space functions; the total angular momentum of the 
centroid dynamics and the Hamiltonian itself. Therefore, the nonlinear dynamics of centroid 
motion is completely intcgrable. We exploited the conservation of the angular momentum to 
reduce the dynamics to that of a 1-degree of freedom system, where the angular momentum 
appears as a constant parameter; specifically, we used the quantity C, which is the amount 
by which the angular momentum is lower than the maximum value which is attained for 
discs with zero centroid eccentricities. The results from a preliminary exploration of the 
nonlinear dynamics are given below. 
I. Case £ = 0: 

(i) A special case includes the initially zero-eccentricity discs, whose linear instability was 
discussed earher. We followed the global phase space structure, as a function of varying mass 
ratio, II, at fixed softening. When /i < /icrit, the zero eccentricity state is, of course, stable. 

(ii) In addition, there are two other equilibria, one stable and the other unstable, both 
of which correspond to eccentric ± discs which processes steadily. The unstable (stable) 
equilibrium corresponds to the discs having aligned (anti-aligned) periapses. The stable 
equilibrium remains stable for all values of /i . 
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(iii) When n exceeds /icrit, the zero-eccentricity equihbrium goes unstable by merging 
with the unstable equilibrium. 

II. Case £ ^ 0: 

(i) We followed the global phase space structure, as a function of varying £, at fixed mass 
ratio and softening. These fixed parameters were chosen such that the zero-eccentricity discs 
would have been in the stable regime; specifically, we chose = r = 0.1 . 

(ii) When £ = 0, there are three equilibria, two stable and one unstable, as discussed 
above. 

(iii) With increasing £, the zero-eccentricity equilibrium remains an equilibrium point, 
but now corresponds to steadily precessing eccentric ± discs with aligned periapses. At the 
critical value of £ ~ 0.0271 , it merges with the unstable equilibrium. 

(iv) Meanwhile, the other stable equilibrium remains stable as C increases. However, near 
C ^ 0.0031 , the steadily precessing eccentric ± discs switch from anti-aligned to aligned 
periapses. Then, near C ^ 0.0433, the equilibrium becomes unstable and, in a pitchfork-like 
bifurcation, there emerge two stable and non-aligned equilibria. The remarkable feature of 
these new equilibria is that they are neither aligned nor anti-aligned. This suggests that, 
for large enough values of £, the stable, uniformly precessing counter-rotating discs have 
periapses that are neither aligned or anti-aligned. 

(v) Increasing C still further maintains the equilibrium structure with the stable equilibria 
increasing in eccentricity and misalignment, in a fashion reminiscent of the displacement of 
an unstable bead on a rapidly spinning hoop. 

The brief report above is but a preliminary account of the vast and rich dynamics of 
this gravitationally coupled system. Here, we have attempted a self-contained presentation, 
whose aim is to point out the variety of steadily precessing eccentric configurations that 
are allowed and how their properties and stability depend on parameters such as the disc 
mass ratios and angular momentum. Straightforward generalizations are possible for discs 
with different values of the semi-major axes, or possibly a range of values of the semi- 
major axes. Razor thin ±-disks, with a spread in semi-major axes, are already known to 
exhibit some of the dynamical properites of the singular distributions studied here: a- a 
linear stability thre shold which r e flects an interplay bet ween softening (heat) and ±-mass 
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201ol): the possibility of a non-aligned. 



ratio (self-gravity) (ITouma 
uniformly precessing equilibrium, flrst identifled by ISambhus fc Sridhad (120021 ) as a promis- 
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ing stellar dynamical model of the double nucleus of M31. A multi-centroid generalization 
of our two-centroid theory, one which allows for a range in the semi-major axes, is expected 
to provide idependent confirmation of these results, and might, in addition, exhibit new and 
unsuspected dependences on the semi-major axes. Other generalizations involve the (com- 
bined) effects of a central density cusp and general relativistic corrections, both of which 
induce apse-precession that can be strong enough to alter the stability threshhold dramati- 
cally. A limitation of the work in this paper is that we have not been able to discus s quest ions 



regarding the saturation of instabilities, such as those explored in iTouma et al.l (120091 1. To 
attempt such a description would involve taking any of the time-dependent DFs of this pa- 
per as an "unperturbed" solution, and solving the linearized CBE for perturbations about 
this state. Such a saturated state could well describe the double-peaked, lopsided nuclear 
disc of M31. 

Whereas we have focused attention on counter-rotating nearly Keplerian discs, our 
general formulatio r i of t he CBE for nearly Keplerian star clusters extends the work of 



Sridhar fc Touma 



(I1999I ) to self-gravitating systems. In particular, our formalism of § 2 
applies to fully three dimensional clusters. In addition to the eccentricity-periapse dynamics 
which applies to planar discs, we will also have to consider the inclination-node degrees of 
freedom. An expansion of the Hamiltonian to fourth order in the (sine of the) inclinations 
will need only a modest extension of the methods introduced in this paper. Then, it will be 
possible to explore the stellar dynamics of counter-rotation, eccentricity and inclination all 
considered together for nearly Keplerian systems. 
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APPENDIX A: RING-RING INTERACTION FUNCTION TO 4™ ORDER 
IN THE ECCENTRICITIES 

We provide exphcit expressions for the expansion of the orbit averaged interaction function 
between two softened, coplanar Gaussian rings, given in equation ([5]), up to 4th order in the 
eccentricities. The expansion is carried out for arbitrary semi-major axes with the help of 
classical techniques of celestial mechanics, generalized to softened interactions. These same 
methods can be used to recover expansions to arbitrary orders in eccentricity and inclination. 
D etails of the techniques, ex pansions, accuracy and conditions for convergence are discussed 



m 



Mroueh fc Toumal (120111 ). 



We consider two coplanar rings with orbital elements {a,e,g) and {a',e',g'), define 

p = min(a, a')/max(a, a') ; r = 6/max(a,a'), (Al) 
and denote the softened analog of the classical Laplace coefficients by, 

2 r cos{mt)dt 
Jo [1 + p2 + r2 - 2p cos(t)] 
The 4th order expansion (with the constant ignored) takes the form: 



B7{p.r) = -I . . . : : (A2) 
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The expansion in the text, equation (!T9|) . is expressed in terms of the eccentricity vectors: 



e = (ecos(7, esingf) ; e = (e cos (7 , e sin^f ), (A12) 

and neglects terms that are independent of e. With these constraints, equation flA3|) reduces 
to: 



-{c°oe' + c}ie ■ e' + c%e^ + c^^e^ (e ■ e') 



2max(a, a' 

+ (4 - 4)e2e'2 + 24 (e ■ e')' + c\, (e ■ e') e'^} . (A13) 

By further speciahzing to rings of equal semi-major axis, i.e. a = a' = ao, p = 1 and 
r = 6/ao, and noting that in that case c\t^ = c\^, we can identify coefficients (Eqn. [19]) with 
coefficients in the expansion above: 
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We end with a brief remark about the convergence of the series expansion. The unsoftened 
expansions in powers of eccentricity are known to have a finite radius of convergence in the 
case of non-intersecting rings, in addition to blowing up when rings intersect. The later is 
alleviated by softening interactions as we have done. However, convergence of the softened 
series remains an issue: in the overlapping configurations considered here {a = a' = oq), and 
for a given eccentricity e, the softening has to be larger than a critical value be = a^e for 
the series to converge. Mo re details on conve r gence analysis and accuracy of the 4th order 



expansion can be found in 



Mroueh fc Toumal (1201 If ). 



